{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's play a bit with our current XHM results. First, let's check how the variants we found vary with data quality. We could assume Q=60 like their paper, but it's worth checking how it changes across quality thresholds. Also, let's plot that value per samples, as CCGO samples seemed to have some noise to them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "import vcf\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "vcf_reader = vcf.Reader(open('/data/NCR_SBRB/simplex/xhmm/DATA.vcf', 'r'))\n",
    "SQs = {}\n",
    "LQs = {}\n",
    "RQs = {}\n",
    "for sample in vcf_reader.samples:\n",
    "    SQs[sample] = []\n",
    "    LQs[sample] = []\n",
    "    RQs[sample] = []\n",
    "for record in vcf_reader:\n",
    "    for sample in record.samples:\n",
    "        SQs[sample.sample].append(sample.data.SQ)\n",
    "        LQs[sample.sample].append(sample.data.LQ)\n",
    "        RQs[sample.sample].append(sample.data.RQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "quality = np.arange(10, 100, 10)\n",
    "cnvs = np.zeros([len(quality), len(vcf_reader.samples)])\n",
    "for q, qual in enumerate(quality):\n",
    "    qual_prime = qual / 2.0\n",
    "    for s, sample in enumerate(vcf_reader.samples):\n",
    "        some_deletion = np.any(np.array(SQs[sample]) >= qual, axis=1)\n",
    "        good_left = np.any(np.array(LQs[sample]) >= qual_prime, axis=1)\n",
    "        good_right = np.any(np.array(RQs[sample]) >= qual_prime, axis=1)\n",
    "        idx = np.logical_and(some_deletion,\n",
    "                             np.logical_or(good_left, good_right))\n",
    "        cnvs[q, s] = np.sum(idx) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2aaafea2c410>"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm4nXV17z/rzCcJQ8IYklBAIjIoUXMpVUvRSEXKldZW\nb7zFYuGKz71UsXqvgm3V1nIfbiuotZY2AlpHShFaLo8iQ+ECDgyBABkIhCQkOSQ5SSDTGfa47h/v\ne8LO8Qzvfve7f2sP65Pnfc4+79n7fH97n52117t+axBVxXEcx2lcOqwX4DiO40yNG2rHcZwGxw21\n4zhOg+OG2nEcp8FxQ+04jtPguKF2HMdpcNxQO47jNDg1GWoROV9E1orIOhG5KqtFOY7jOK8haQte\nRKQTeB44D9gCPA58SFVXZ7c8x3Ecp6uGx54FrFPV9QAicgtwETCpod73JxcEL4Pc/BMJLQnAw7nZ\nJrrPdhWCaz5X2h1cE2BL7hUT3X2FYRvd/IiJbqFcMtEdGXmp5v+8hZ3rE9mc7iNPsjEUCakl9DEP\n2Fzx/Zb43EGIyOUi8oSIPPGtVZtqkHMcx2lPavGoE6Gqy4BlADvO+y0dXTtUb8mDeHHouKB6Ywz0\n2fRQ2VjeH1zTyrPdldtrojtcyJnoFo0826buB2T0mmVNLYZ6AFhQ8f38+JzjOE5jUCparyATajHU\njwMLReREIgO9FPivUz1g37beGuTSsbejM7gmwG7Cx4oBCloOrvlrvUeyYXRHcN0usfnb9nTU/UJ0\nQgSbMGq+ZPNezgI1+P9QD1K/41S1KCJ/AvwU6ARuVtVVma3MaRosjLTjJKLc5oYaQFV/DPw4o7U4\njuNkS7t71Gk49vKTQ8oB8Ja/2xhcE2Df6BwT3T2dM4Jr7u85NLgmgGKzyZUzCgV0dhgVEjfzfpxv\nJjqO4zQ47lGDiGwE9hF95hZVdfGU9z8xvEcdLTE8I0bp80Ma3tvbkbdJk9s5ssdEd9TIoy63iNEJ\niXrWxwHeqao7M/g9juM42eKbidWjG9aFlANgNGcT3bFqS3iU9AXXPKbnsOCaALly3kS3p2iju78w\naqJrVWiTCRlehYjIzcCFwKCqnhGfmwP8C3ACsBH4oKq+Gv/sauAyoojDJ1T1p/H5twLfBvqJkjGu\n1Gmqimq1JwrcJyLLReTyGn+X4zhOtpRLyY5kfBs4f9y5q4D7VXUhcH/8PSJyGlFtyenxY/4hbmQH\ncAPwUWBhfIz/nb9Cre7mO1R1QESOBu4VkedU9aHKO8QG/HKAr559Ch95/a+0A6krQ0bNkV7us/FC\nNpfDlugDHNbZx+b8q8F1h4w8zBEjj7pkdBnf1CXkGXrUqvqQiJww7vRFwLnx7X8GHgQ+G5+/RVVz\nwAYRWQecFe/rHaqqvwQQke8Avwv8ZCrtmjxqVR2Ivw4CdxB11Bt/n2WqulhVF4c20k4YLIy04ySi\nVEx2pOcYVd0a394GHBPfnqxp3bz49vjzU5LaoxaRmUCHqu6Lb/828FdTPab7N85IK5eaY5c/F1wT\nYMGeo0x0RzpnBdc8qW8WKwrh95N3d4ePxwOUjTzMUoeNR93UMeqEVyGVV/4xy+KGcolRVRWRurw5\nagl9HAPcISJjv+cHqnp3JqtymgoLI+04SVBN9iFT2eWzSraLyFxV3Soic4HB+PxkTesG4tvjz09J\n6tCHqq5X1TPj43RVvSbt73Icx6kLWk52pOdO4JL49iXAv1ecXyoivXHjuoXAY3GYZK+InC2Rl/tH\nFY+ZlLC5awYJ+51dNpeLM43SN3Od4YX7O3qCawL0GunmO22KKPbmbSbLWG1iZkKGaxeRHxJtHB4p\nIluALwDXAreKyGXAS8AHAVR1lYjcSjTxqghcoa+59/+D19LzfsI0G4ngJeSO47Qy2WZ9fGiSHy2Z\n5P7XAL8SaVDVJ4CqNuymNdTVJnlP+btOWFjN2jLhkHnPBNcEOHmbTerYuu7wXuag2Hzed4pNWVGH\n2PQHsHq+JZrYo27iXtqVJPnLf5uESd6O4zgNRbmc7GhwpnWFqkzynpodW6e9S9YMPj8zuCbAI702\nqWMbNHzBy7ZS+DmNYOfZWrU5zRs1GLJqJ5sJLdLIKu0162RJ3o7jOI1DE3jLSag5uDhdkndlIvnX\nP7yEy855U62SVWGV9dFjM6iagoR/vvO7DmV7MbxXbVV40mU0h3NGd/iZo2A3dT0T2txQT5bk/StU\nJpKP3PipJr6GcibDwkg7ThK0RTYT0xrqsSTvazk4yXtqsf/8sZRy6Zn38pRV7XXjXTfuM9HdKeFj\n8id3z2RFKXy/j0Gj2KnVNHCrmHxT0y4x6mqSvJ32xMJIO04i2iX0UW2St+M4TsPQLh51luT/KXwY\n4tl/HAmuCfBQ9yEmujsIX2iTT9j4JmtyJZu+0IWyTZrcaNEm3irNHHJpF4/acRynaWkXj3qSEvIv\nEo2S2RHf7XOq+uNpf9eM/vQrTcmh/buDawJ0FGaY6L5CeC/z1ZJNsyArCkb9mUtGRqfUzP2oi60x\nhTxtCTnAV1R1UXxMa6Qdx3GCU/82p0FIW0KeUi18oUBvv80naqdR+qZF2x4rT69slJ5n1RypbPU6\nN/PMxBaJUdfyjvu4iDwjIjeLiM0EWcdxnKloF496Em4AvgRo/PU64NKJ7lhZQn7dGQu55Pi5KSXT\n8eKOsHpjPN9vk5HwUj58TL6no4stw+HHcVllX5g1R2pmz9aKFvGoUxlqVd0+dltEvgncNcV9D5SQ\n7/qd3/J3WgtiYaQdJxFN4C0nIZWhHuvzEX/7e8DK7JbkOI6TES2S9ZG2hPxcEVlEFPrYCCRq4tEx\nI3zivBhtOB1ulKI+uzN8WmC+7/DgmgC78ntNdPfnbab3OClokXBR2hLym+qwFsdxnGxp5xh1Wmbd\ncHNIOQB+47r/FVwTYMGPbLrn9eeOCK75s26bEuO80WYiNrVM7M3ZFBZZbZ5mghtqx3GcBqddNhNF\nZAHwHaJxWwosU9WvpZlEXnzollrXWzWv3GtTQr5y37Emujt7w3s/nUZNe6wKbdptZmJTl5CXmnjt\nFSQpeCkCn1bV04CzgStE5DR8ErnjOI1OG00h3wpsjW/vE5E1wDxSTCL/+X9/uoalpuOH/TZFk6s6\nXjbRHdhrk9O8azR8TL5k9B/Maip32ez5NjFNYISTUFWMOu758WbgUXwSuRNjYaQdJxHtEqMeQ0Rm\nAT8CPqmqeyubiU81ibyyhPzKQxZzQf/raltxlVj9mbpN2iNBV0f4/eFjZsxmuBg+t9gqn7nYzDHb\nNkPLTX09cIBE1kREuomM9PdV9fb49PZ4AjlTTSJX1WWqulhVF4c20k4YLIy04ySiRWLU0xpqiVzn\nm4A1qnp9xY/GJpFDFZPIHcdxglEqJTsanCTXyW8HPgw8KyIr4nOfI8Uk8rcnKjTPlvI/2aSKP9Zn\nk573SEd3cM3dJZu5lJsPDBgKy1DB5gqiYJQGadUHOxOawFtOQpKsj0eAyd4hPonccZzGpV0MdZas\n/VZ4z2tVr8008OexKffdWdwfXHMwZ1NUZOXptVvhSVNvx7VIUyab1ATHcZwQZLiZKCJ/KiKrRGSl\niPxQRPpEZI6I3CsiL8RfZ1fc/2oRWScia0XkPbU8jVpKyL9IlZPIT/2bRbWsNRVHXf9IcE2AYwZs\nJsvM7A2vu7qjN7gmwPrhbSa6VgUvYhSjburJMhml54nIPOATwGmqOiIitwJLgdOIKrSvFZGriCq0\nPxtXby8FTgeOA+4TkderaqrLoiShj7ES8idF5BBguYjcG//sK6r65TTCjuM4dSfbjI4uoF9ECkQ9\nFF8GrmbiCu2LgFtUNQdsEJF1wFnAL9IKT8kUJeRVU165Os3DamLn9lnBNQGGjIJKnZPu+9aPN3Yf\nwVOF8BkYVm1OZ3TZXEEUjGLUVrMps0Az2kxU1QER+TKwCRgB7lHVe0RksgrtecAvK37FFlLaTagy\nRj2uhBx8ErkDJkbacRJR1kSHiFwuIk9UHJdX/prYvl0EnEgUypgpIhdX3kejGFFd4kSJDfX4EnKi\nSeQnAYuIPO7rJnncgRfg5qfWZ7Bkx3GchGg50VFZQR0fy8b9pncDG1R1h6oWgNuBtzF5hfYAsKDi\n8fPjc6lIlJ43UQl50knklVPIh//20uC7Ekcc/VJoSQBK223SAvdp+MtUq0vjWV19Jrp7yjapl1b9\nt6269mVCdr0+NgFni8gMotDHEuAJYIioMvtaDq7QvhP4gYhcT+SBLwQeSyueJOtjwhJyn0TuOE7D\nU8wmrq+qj4rIbcCTRAkWTxE5oLOYoEJbVVfFmSGr4/tfkTbjA2orIf9Q1ZPIZ81Mt8oaGNnfE1wT\noGCTSUXeoF9g2ShdzcrDtEpX6zBKz+voaOJyiwzfI6r6BeAL407nmKRCW1WvAa7JQruWEvIpc6Yd\nx3HMaZE2p0FLyMsvbgopB8Duvf3BNQH2d9u8QQZLQ8E1c2WbGYJWHnVXZ6eJrtnU9SYmq/Q8a3wK\nueM4rUu7eNQi0gc8BPTG979NVb+QZgp5aVd4b29/6YjgmgBbemy8TAuO7j6ULblXgut2ik3stGjU\nv9hsRqSXkJuT5J2eA96lqmcS5UyfLyJn41PInRgLI+04iWiXwQFxtc1Y78zu+FBSTCHv+/O/TrnM\n9Jw9//rp71QHjvquTTP9xy1mDPcew4Od4XOLV+UnnP5Wd0ZKeRPdHrWJVFqVrmdBu81M7IxT8waB\ne1XVp5A7B7Aw0o6TiIQl5I1OIkOtqiVVXURUBnmWiJwx7ueT1rhXlpDfeMu/1bxgx3GcxLTIcNuq\nrqVUdbeIPACcT1zjrqpbp5tCTlxCPvKtz2j5odsnulvdGH1y+/R3qgMbRlM3yqpNty/8m25v2SgU\nIDahgCN7DzPRHeq0CaftGt1nopsJTeAtJyHJFPKjROTw+HY/cB7wHD6F3HGcRqdFQh9JXJK5wD+L\nSCeRYb9VVe8SkV9Q5RRyjj+5lrWmojT66PR3qgMbe2xSx17Q8DMTN+WnzMqsG1tHbbJNRoo2VxDD\neZvp541vxiZHS40f1khCkqyPZ4h6UI8/vwufQu44TiPTBN5yEsKWkD/8QEg5AJ5eeWxwTYD9fUaN\neww0rZoy9XZ2m+haTSHv6jRKzzN6vlnQKul5XkLuOE7r0i6GeooS8i9S5RTyjpNPqm21KTj1hFSz\nJGtm9dbjTHT7OsI3DHpL31xeLoXPpbZqUvRqOfw+AEDZqAlVU9MiL1kSj3qshHx/POnlERH5Sfwz\nn0LumBhpx0mCFlvDUtdSQu44jtPYtIadTjwzsRNYDpwMfCMeS/Neoinkf0Q0O+zT03XP0z17al1v\n1bw6GH6qDMChRu0RhiW88J6STSHG3kL4boxg17XPbPMU30y0ppYS8uqnkD+yKqNlO47jJKCc8Ghw\nUpeQV8amk04hf/GM9+jWZ1+oYbnV87PiUUH1xljenTPR3VkMHy8ezIe/UgIoG/VJtuomlzMqtGlm\nn7RtPOrJSsjj/h5j+BRyx3EajzbyqCcrIf9utVPIj353+Ingb/s3G29vdu5QE92f984Jrtk/wyZ2\nuilvU0JuNQ083x3+/w/ASMHGk88Cbd7w+kHUUkL+4bqsyHEcJyNaJfU8aGXiyLN7Q8oB8PPSguCa\nAE/12sSoXyyFf40BXhjeOv2dMmbXiE37zaJRjFqNosXNPTPRegHZ4CXkTs1YGGnHSYJ71CnonBFS\nLeKEvE2Qan2/zWdgv4SPF79p5vG8lN8VXLcDm1jxcNHmasmqvWpzz0y0XkE2JM7cj+cmPiUid8Xf\nzxGRe0Xkhfjr7Pot02lkLIy04yRBS5LoaHSqKbG6ElhT8f1VwP2quhC4P/7ecRynYdBysqPRSVpC\nPh/4HeAa4FPx6YuAc+Pb/ww8CHx2qt/T96aj06yxJg553Cb0saDUZ6K7waB73iGd/cE1AfYYlZBb\n9aMuGVmUZt5M1HLje8tJSOpRfxX4DAfvoR6jqmO7SNuAYyZ64EEl5E+tT79Sx3GcKmkbj1pELgQG\nVXW5iJw70X1UVUVkwo/dyhLyoT/7QPCP5s4Om79CzuiDvI/wHrWVp2dFV2f41xgMX2ebHlSZoNoa\nHnWS0MfbgfeJyAVAH3CoiHwP2C4ic1V1a1xOPljPhTqO41RLq/gQSSoTrwauBog96v+pqheLyN8C\nlwDXxl//fbrfJfPDzy/s6tgcXBOgaPRBvl8LwTVHykZTuY3S5IYLNrr5Yvi/bbNTzjCjI+55dCNw\nBlHrjEuBtcC/ACcQtdL44Fi7ZxG5GrgMKAGfUNWfptWu5aLmWuA8EXkBeHf8veM4TsOgZUl0JORr\nwN2q+gbgTKIsuAmz30TkNGApcDpwPvAPcb+kVFTb5vRBouwOVHUXsKSaxxceWzP9nTJmf8EmvXu4\ny2anXAyKQF7XeyQvF8KXrs/qssmsKZZsCkAKYpNt4lkfICKHAecAHwFQ1TyQF5HJst8uAm5R1Ryw\nQUTWAWcBqYa4NvE2gdMoWBhpx0mCarKjMjstPi4f96tOJBrk/a248O9GEZnJ5Nlv84DKuOuW+Fwq\ngtY59/z+BSHlADjx6SkHo9eNU7dNmK1Ydzb3hC8hP6znCFaMbguuO1KyiY1bZV+IUXvVZiapR12Z\nnTYJXcBbgI/Howi/xrgiv6my32qllhLyL4rIgIisiI/wVthpCCyMtOMkQVUSHQnYAmxR1Ufj728j\nMtzbx4aojMt+GwAqW3fOj8+lopYScoCvqOqi+LBxXR3HcSahVJJEx3So6jZgs4icEp9aAqwG7iTK\neoODs9/uBJaKSK+InAgsBB5L+zxqKSGvmuF/mjaDL3PWDMyd/k51YL3NPhc7y+FTx/JGYzQsNk7B\nbgp5h5FuuYmbOmdc8PJx4Psi0gOsB/6YeOqViFwGvAR8MNLVVSJyK5ExLwJXqGrqXeikMeqxEvJD\nxi9cRP4IeAL49Fj+YCVxUP5ygOveuJBLfs3GcDqO035k2etDVVcAiyf40YTZb6p6DZFzWzO1lJDf\nAHyJKPH7S8B1RAng4xd7IEi/7ZxztRB4hOFAt01f6A0yaqK7rxRetyt9emhNWE08KbfZpJVSuZk9\nausVZEPqEnJVvXjsDiLyTeCuOq3RcRwnFa3SPa+WEvK5FfmDvwesnO53zVw0q4alpuPYtTbFCYd2\n2XjyvQYTXorpQ2810WXQ0hWgp9Pmb2vl2VrNiMyCUrk1SkVqecf9jYgsIgp9bAQ+lsmKHMdxMqKd\nQh8HGFdC/uFqxXJrw0+NLsihwTUBDlebT/LDO3rDa/bN5fncjuC6vR09wTUBCkaDAzo7bN5TVlcu\nWVBuozanjjMlFkbacZLQTv2oHcdxmpK2Cn2IyEZgH1Ff1aKqLhaROUzSh3UyZlxwai1rTcUbNmwK\nrgnw0ugcE90+g8vUGZ3hwy0Ae4sjJrojRZseI2azGpt4M7FVQh/VBL3eGZeKjyV8+xRyx3EamlK5\nI9HR6NQS+qh6CnnH286rQS4dx/7874NrApzyoI33s6cvfO36LqMp5CPdNp5t2Sgd8ZXcfhPdfIZT\nUkLTIpGPxB61AveJyPKKPq1VTyG/6Y57a1yu4zhOcsoqiY5GJ6lH/Q5VHRCRo4F7ReS5yh8mnUI+\n/NWPafln99S04GpZ97PDguqNcV+/zR//F4Wt098pY3YXh4JrAgyO7jbRtYsV2xS8lJt4QmyrZH0k\n8qhVdSD+OgjcQTRSZrI+rI7jOA1BOeHR6CRpyjQT6FDVffHt3wb+itf6sCaeQr71xo01LTYND3Yc\nFVwT4IWyjbeXK4efVN3f0cOOfOBuW9g1KRo1yvpo5uZIVqhRK9ysSRL6OAa4Ix4D1AX8QFXvFpHH\nmaAPq9N+WBhpx0lCsUVCH0maMq0nGo0+/nzVU8h/+Up47/bZHptc21GjSdWjBh71IV0z2GEQLx4t\nhn+uYFdSbTUooalj1G3kUTvOlFgYacdJQvN+xByMG2rHcVqWtvKoJykh/yLwUWCsI8/nphtw+/5v\nnJF+pSk5+y8eCK4J8LNXbTYxV/SHT0dc2WPjUW/K7TLR3ZO3SUe02sQcLdmEmLKgHT3qd6rqznHn\nvqKqX85yQY7jOFlRaiePOiu2/+X9IeUAuHfPhAWTdeehLhuva31+yr5YdWF7LrwmwLZhG0++mZsU\npcEqDTILWmQSV00l5BBNIX9GRG4Wkdl1WJ/jOE5qykiio9GppYQ80RTy2LBfDvD5I97IBw49PpOF\nJ2W/0d9gSG3ienG+e1BGSjaxU6vZhSNGHnW5iT1bK1rlFUtdQq6q21W1pKpl4JtEZeUTPXaZqi5W\n1cWhjbTjOO1N25eQp5lCfvyS8J7XGXfaNNDZ3B9+4jrA3nIuuOaps+YzWNgbXHe32LT9tMJqGnje\nqLAoC8oGV5j1oJYS8u/6FHIHMDHSjpOEVtn2raWEvOop5I7jOCFplayPoLsxHUceElIOsMuj7DZR\ntaFgdEluFgow6kdtlSbXzBtyzZDRkQQvIXccp2Vp5g+ZSpKWkB8O3AicQfTcLwXWUuUUcn01fBHI\ncEf4GYIAOxg10d1TCt8tsGDQsS/SbZUIZDKsriCamVYJfSQtePkacLeqvoEoXr0Gn0LuOE6D007p\neYcB5wAfAVDVPJAXkaqnkN9615HpV5qSX/TaeLaDpWETXYs4plXTng6j1CurWHG7Pd8syHqAuoh0\nAk8AA6p6oYjMYZLIgohcDVxGlHzyCVX9aVrdJB71iUQd8r4lIk+JyI1xPnWiKeSO4zhW1MGjvpIo\nojDGhJEFETkNWAqcDpwP/ENs5FORJEbdBbwF+LiqPioiX2NcmGOqKeSVJeR/MeeN/MEhv5Z2ranY\noTZTyMs2Q0AoGFzIHdE9iw0jO6a/Y8bsLdg0vrIo0wcol5vXs7Uiy/8NIjIf+B3gGuBT8enJIgsX\nAbeoag7YICLriKq3f5FGO4lHvQXYoqqPxt/fRmS4E00hrywhD22knTBYGGnHSYJKsiMhXwU+w8H2\nf7LIwjxgc8X9tsTnUpGk4GWbiGwWkVNUdS3RnMTV8VHVFPLjfz18ye+bHpwZXBPg5X4bl7o78f5w\ndry+/xheHA1vrHs7bbLVcx1GWS5G+dvNTFKPuvLKP2aZqi6r+PmFwKCqLheRcyf6HVNFFmolaR71\nx4Hvi0gPsB74YyJv3KeQOyZG2nGSkDShMTbKy6a4y9uB94nIBUAfcKiIfI84sqCqW8dFFgaABRWP\nnx+fS0XS7nkr4vDFm1T1d1X1VVXdpapLVHWhqr5bVV9JuwjHcZx6UJZkx3So6tWqOl9VTyDaJPwP\nVb0YuJMoogAHRxbuBJaKSK+InAgsBB5L+zyCVibuWx/+svzFHpviy83lPTa6BhNecmWbftT78uGL\newBKapN529VhE05r5ok2Af5S1zJBZEFVV4nIrUQh4iJwhaqmfiG9hNxxnJalHoZaVR8kyu5AVXcR\n7dtNdL9riDJEaqaWEvL3UOUU8mI+vEe9t8MopclItq8j/AZbzkAToKNFGu4kpZk9WytaJaExqUc9\nVkL+B/GG4gwiQ+1TyB3HaVhapddHLSXkVYsd+6m3Vv2YWvn9a54Jrgkwe2SOie6Tvb3BNZ+WCVPo\n606xbJOutr9o05bAKh1xuBB+alBWtMo1SC0l5OBTyB3HaWDKaKKj0amlhPzvqXIK+df/cAmXnvPG\njJaeDDGKUeeMLrlGDTISTuk+kg3F3cF1Z3b1B9cEEAm/1wIwVLDJcrEqmc+CZuiMl4TUJeRpppCH\nNtJOGCyMtOMkQRMejc60hlpVtwGbReSU+NQSYPVYn4+YRFPIHcdxQtI2/ahjJioh/7uqp5Ab9Coo\n5G2KBPbUp+R/WizibYX0efw1MVS0CQUMF20214aMNvWauR910ej/YdYkMtSqugJYPO60TyF3HKeh\naQ0zHbgyUY5/XUg5AA4/alVwTYCTthxuoru9O3wK1xGdM4JrAox225SuW2UJWG3qjRRtXucsaIaw\nRhK8hNxxnJalGVLvkpCk4OUUoplgY5wEfB74DlVOIS898nDadaYmP2zzWbTHJoOLHRq+GONVg8nn\nYFd4YtUMymoKeTNPe28NM50s62Otqi5S1UXAW4Fh4A58CrnjOA1Ou2V9jLEEeFFVX0ozhTy/eme1\n66uZ5TvnB9cEeKHXZgqIxVi947oPZUNuV3jhNmPUKFbczFkfpRbxqas11EuBH8a3fQq5A+BG2mlY\nmsFbTkJiQx3nUL8PuHr8z5JOIf/60nO59B2np1xqOt7wzEtB9cbYn7eZfr6vO3yu+qK+Y9lY3Btc\n95XCvuCaAJ0dNhsQfV09Jrr5Jp7VqC3iUVfzjnsv8KSqbo+/r3oKeWgj7YTBwkg7ThJaJUZdjaH+\nEK+FPWDyWWGO4zgNQTt1zyNua3oeB5eJTzgrbCpKazemWGJtdBp1z3vFpnKdrqo+e7Nhf7m9Spu7\nxOaPWxSbNLlm7p7X+CY4GUlLyIeAI8adm3RWmOM4TiNQbBFTHbQaRIfCe17PDx0VXBNgX6/NG2Rb\nOXwxRp/YTB6x8vR254ZMdK2uIKwKbbKgVTYTvYTccZyWpRk2CpNQSwn54VQ5hbwwGD5hP2fkdW0V\nm+KE4VJ43W15m8EBQwWbEvJOowkvJSOz09HUMeo28ahVdS2wCEBEOoEBohLyP8ankDuO08C0jUc9\njsoS8qrFZnz6kunvlDHn/Pl3g2sClNcfZ6Lb2xN+xvDpXbN5rrQnuK5V9sW+4rCJbq5k05Zgf97m\nyiULSk1c/l5JtddwlSXk4FPIHTAx0o6ThLbKo4YJS8hvoMop5Ne/ZSEfOSmsp9k3Fx54OLx3u7nX\nJq63vrw/uGaPdLJ6+OXgugCDw+Hj4zmr5kgmqs1N28SoKziohLyilBwR+SZw10QPUtVlwDKAVz9w\nbvBXzcJItxvtZKSd5qJVYtSpS8h9CrnjOI1OW4U+Jikh/5tqp5A/8vDc6e6SOb/ss/kjrCnbeHsb\nDVqO7ivYbK51d9hsJhaMuuc1c19oK9oq9DFJCblPIXccp6FplayPoJWJ77wkfAn5ybfbpBbdk5tj\novt0f18zrZMbAAANCElEQVRwzfXds4JrAmzL22SbWF1BWJVyDxVsmm5lQTOENZJgNILVcRyn/mTV\nj1pEFojIAyKyWkRWiciV8fk5InKviLwQf51d8ZirRWSdiKwVkffU8jySxqj/FPhvRPHoZ4mqEmdQ\n5RTyvQ9P+eO68Mx+m6yPlb02XsiGYngvcyAX/u8KsGvUZmDBsJGHWTLyqJvZJ80wRl0EPq2qT4rI\nIcByEbkX+AjRkO9rReQqoiHfnxWR04jqTk4HjgPuE5HXq2qqP+K0HrWIzAM+ASxW1TOAzngBPoXc\ncZyGJqusD1XdqqpPxrf3AWuAecBFRMO9ib/+bnz7IuAWVc2p6gZgHXBW2ueRNEbdBfSLSIHIk36Z\nqPDl3IoFPsg0U8hnLgwfaTnyRZt5b8erzXy7LRK+IeLJfUezrWDg3YYPxwPQITazGgtGHnWp3LzZ\nyEkzZSoL82KWxTUgE933BODNwKNMPuR7HvDLiodtic+lIklTpgER+TKwCRgB7lHVe0TEp5A7ADZG\n2nESUEoY+qgszJsKEZkF/Aj4pKrurex5NNWQ71pJEvqYTeTGn0gUa5kpIhdX3kejj61Jp5CLyBMi\n8sS3nx/IYMmO4zjJyLLgRUS6iYz091X19vj0ZEO+B4AFFQ+fH59LRZLr5HcDG1R1R7yY24G3jS1Q\nVbdON4Wc+JPqyQUX6Ytb0i41HSt6baaPPI9NCteIhu+wtrcYfqoMwJ58+L4mAKNFmy52JbUJQRRL\nNuHDLMiqSEgi1/kmYI2qXl/xo7Eh39dy8JDvO4EfiMj1RA7uQuCxtPpJDPUm4GwRmUEU+lgCPAEM\nTbJAx3GchiDDPOq3Ax8GnhWRFfG5zzHJkG9VXSUitwKriTJGrkib8QHJYtSPishtwJOx4FNEHvKs\niRY4Fad8NHxhRP4fbTzbUvcME93RzvAbTvu7bdLVhos2xUzlTpuENau0wGYmq/Q8VX0EmKwl5oRD\nvlX1GuCaLPSTlpB/AfjCuNM5fAq54zgNjJeQp6C8M3wxRll7g2sCGDldFA0aO3Z32MxIPrzHpnR9\nZ86mdL3LqAlVM08hb5UScp9C7jhOy9JWhnqSEvKrqHIKec+Vf55+pSk5/dU/C64JMHSPzWfg1v7w\nXubrOmbxRHHH9HdsEYolGw/TyrNt6ink7RL6qCghP01VR+KdzKXxj30KudNWRtppLtrKo2biEvIT\nqhV74by/rvYhNfNQIfywAoBH+21yi7cUw5c3d9LBK8Wh4Lqv5GwqIvNlm7ziZo4VW9EqgwOmrUxU\n1QFgrIR8K7BHVe+Jf+xTyB0TI+04SShpOdHR6NRSQn4DcBKwiMiAXzfJ4w+UkP/r3k2ZLdxxHGc6\nVDXR0eikLiFX1e+N3SHpFPLdf/guJXBp9YL7DwmqN8ZznTabiYEr9AHIq00oQIw2uTrF5200C60S\no07yjjtQQh7Xuy8B1vgUcsdxGh1N+K/RqaWE/MZqp5Bv/vnMmhabhs3dRkUCGM23K+eDa1rNENyf\ntykht4ppWiXJNb4Zm5xyE4Q1klBLCblPIXccp6FpBm85CUEDqSd//tSQcgAc+vUng2sC9A4ea6Lb\n1xNed2WnzaiVl7t3m+gODO000R028uSbYbNtMpohoyMJXkLuOE7L0lahj3g0+keJwmTfVNWvisgc\nqpxCvuMby2tabBoe3mZT8LK8J3ysGODZwi4T3e358I2Kto/YTD+3KjwpN/HsQitaJfSRJI/6DCIj\nfRZwJnChiJyMTyF3YiyMtOMkoaya6Gh0kqTnnQo8qqrDqloE/h/wfiYfk+44jtMQtE16HlF+9DUi\ncgTRKK4LiEZxVT2FvMOgSXPBKKdp1KAvNEDOoPjEqueG1SaXlQdmVeDTzJTST79qKJL0+lgD/B/g\nHuBuYAUcnCScdAr59wdfrn3FjuM4CWmnEnJU9SaiCbyIyP8mqlSuegr5/s++P/grcsb3bBoGbdPw\nxT0Au7rCl8znZthM5X41H75TIMCuURvdspFH3cybmO1UQo6IHB1/PZ4oPv0DXhuTDj6F3HGcBqSt\nPGrgR3GMukA09ny3iEw4Jn0qihvDp1MNFm0KT/Z323ghnck+e7PVNGpSNFy0mcpdMvIwrQxKM8fG\nmyGjIwlJQx+/OcG5XfgUcsdxGphmyOhIQtDKxDsfmRdSDoBf9NsUnrxYssmE2JSzKXgZLto0SLKg\nt6vbRDdXtNkLKDdxGbaXkDtOTDsZaae5aIb4cxKCGur3/5fwXuZ/+r82Te0fKB1povv4DJtskzX5\n8J789pxNCflIyeYqzcroFJp4VmOrxKiTZn1cKSIrRWSViHwyPvdFERkQkRXxcUF9l+o0KhZG2nGS\n0DZZH+N6feSBu0VkbOzWV1T1y3Vcn+M4TmpaJY86SejjQK8PABEZ6/VRNaOrwjfv+eXwguCaAI92\nj5joPpvbHlxzV95m4zRXstlcGy7YpAVahSCaeTOxGbzlJCQJfawEflNEjhCRGUS9Psas38dF5BkR\nuTmeVv4rVJaQf2fL1onu4jiOUxdKWk50NDpJZiauEZGxXh9DvNbr4wbgS0Q9Pr4EXAdcOsHjD5SQ\nbz/3XC0FdkbKRrn6M7GZ1djf0RNcc1ZXf3BNsPOorQp8cmrzfJuZttpMVNWbVPWtqnoO8CrwvKpu\nV9WSqpaBbxLFsB3HcRqGLDcTReR8EVkrIutEJGj//aQTXo5W1cGKXh9njzVkiu/ye0QhkinJD4X3\nMq0uarqMZkbP7gzv3e4u2jS+EqPXuN1K15uZrCoTRaQT+AZwHlFTusdF5E5VXZ2JwDTU0uvj6yKy\niCj0sRH4WJ3W6DiOk4oMNxPPAtap6noAEbmFaHhK4xjqSXp9fDj75TiO42RHhjHqecDmiu+3AL+e\n1S+fjqCViQsevz/1taqIXB5vTFbFR9MK1qjbbJqu27qa7ag7RjE/kMjmiMjlwOUVp5ZZrns8NtvX\n6bh8+ru0jG47Pdd2022n52qpWxWqukxVF1cc4430AK+lJQPMj88FoZkMteM4jhWPAwtF5EQR6QGW\nEg1PCYJ3z3Mcx5kGVS2KyJ8APwU6gZtVdVUo/WYy1FbxIgvddnqu7abbTs/VUjdzVPXHwI8ttKVV\nauEdx3FaFY9RO47jNDgNaajjJk+DIrKy4twcEblXRF6Iv07YBKoGzQUi8oCIrI77bl8ZSLdPRB4T\nkadj3b8MoRtrdIrIU2NtawNpbhSRZ+Me5k8E1D1cRG4TkedEZI2I/EaAv+0pFf3aV4jIXhH5ZADd\nP43fSytF5IfxeyzEazxR3/q667YDDWmogW8D5487dxVwv6ouBO6Pv8+SIvBpVT0NOBu4QkROC6Cb\nA96lqmcCi4DzReTsALoAVwJrKr4PoQnwTlVdpKqLA+p+DbhbVd8AnEn0vOuqq6pr4+e5CHgrMAzc\nUU9dEZkHfAJYrKpnEG18La2nZqxb2bf+TOBCETm53rptQ9KmJaEP4ARgZcX3a4G58e25wNo66/87\nUV1/MF1gBvAkUcVTXXWJ8kDvB94F3BXqNSZqN3DkuHP1fq6HARuI92Qs3lPAbwM/q7cur1XQzSFK\nFrgr1q73a/wB4KaK7/8C+Ezo/7etejSqRz0Rx+hrTaC2AcfUS0hETgDeDDwaQjcOQawABoF7VTWE\n7leJ/iNVdvoJ8RorcJ+ILI+rwULongjsAL4Vh3puFJGZAXQrWQr8ML5dN11VHQC+DGwCtgJ7VPWe\nemrGTNa3PuRr3LI0k6E+gEYfz3VJVxGRWcCPgE+q6kGjS+qlq1G72EVEXu5Z8WVk3XRF5EJgUFWX\nT7Gmer3G74if63uJwkvnBNDtAt4C3KCqbybqq37QJXid31M9wPuAfx3/szr8bWcTNQs6ETgOmCki\nF9dTM/6da4CxvvV381rf+rrqtgvNZKi3i8hcgPjrYNYCItJNZKS/r6q3h9IdQ1V3Aw8Qxefrqft2\n4H0ishG4BXiXiHyvzprAAY8PVR0kiteeFUB3C7AlvlIBuI3IcIf6274XeFJVx+ak1VP33cAGVd2h\nqgXgduBtddYEJu5bH0K3HWgmQ30ncEl8+xKiGHJmiIgANwFrVPX6gLpHicjh8e1+orj4c/XUVdWr\nVXW+qp5AdEn+H6p6cT01AURkpogcMnabKHa6st66qroN2Cwip8SnlhC1p6yrbgUf4rWwB3XW3UTU\nL35G/J5eQrRxWvfnKiJHx1/H+tb/IIRuW2AdJJ/oIHpTbyXqf70FuAw4gmjz6wXgPmBOxprvILos\ne4bosm0FUZyt3rpvAp6KdVcCn4/P11W3Qv9cXttMrPdzPQl4Oj5WAX8W6rkSZdQ8Eb/O/wbMDqQ7\nE9gFHFZxrt6v818SfdivBL4L9AZ6rg8TfQA+DSwJ+T5u9cMrEx3HcRqcZgp9OI7jtCVuqB3HcRoc\nN9SO4zgNjhtqx3GcBscNteM4ToPjhtpxHKfBcUPtOI7T4LihdhzHaXD+P6dQvg6MZy4zAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2aaafed07650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(cnvs.T, xticklabels=quality)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There doesn't seem to be many outliars, which is good, but Q=60 seems a bit optimistic. We'll need to play with it a bit. Not sure if the control subjects are biasing this... need to look at the papers to see if it's better to run within groups? Before we do that, let's restrict the colors first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2aab0a467350>"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD8CAYAAACihcXDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUXGWV9/HvTjqhkwCBAIncNGFk4I2MgCLDADoIyCAy\ngL4OC2ZAVCSMIoI6CwGXiAIOOMiI3JwQLkEQDMjtjcgtioCOkQQC5kIEQxJyISHEXMi9q/f7xzkN\nRZtOn+46Vfuc6t+HdVZ3Vxf1Ox3C00895+xnm7sjIiKx+kWfgIiIaDAWESkEDcYiIgWgwVhEpAA0\nGIuIFIAGYxGRAtBgLCKSgZndbGZLzWx61WPDzOwxM3sp/bh91fcuMLOXzWy2mf1Td69f02BsZken\nQS+b2fm1vJaISMHdChzd6bHzgUnuvicwKf0aMxsNnAS8L/13rjez/lt68V4PxukLXwd8HBgNnJye\ngIhI03H3J4HlnR4+Hhiffj4eOKHq8bvcfYO7vwK8DBy4pddvqeHcDgRedvc5AGZ2V3oCM7v6FzYt\nm6NyPym9jT+6MCT3lltq+d+19yb3WxuS+5N591qtr5F1zBm409+cCYypemisu4/N8K+OcPfF6eev\nASPSz3cFfl/1vAXpY12qZZliV+DV7sLMbIyZTTGzKeNuu7OGOBGR+nD3se5+QNWRZSDu/BoO9HrC\nWfdftekPNRY0M5Z8zfvIF0Nyr1o7NCT3ufaFIbmrNsbMjHPRXql3whIz29ndF5vZzsDS9PGFwO5V\nz9stfaxLtcyMexwmItJQlbZsR+89CJyWfn4a8EDV4yeZ2VZmNgrYE/jDll6olpnxM8CeadBCkiuH\n/1rD64n0yIiP1Lyk2CtH3R+zdrt9684huXdv/FNIbh7c23N7LTO7EzgM2NHMFgDfBi4HJpjZ6cA8\n4MQk12eY2QSSa2htwFnuvsVpeq//Vrl7m5l9GXgE6A/c7O4zevt6IiK5a89vMHb3k7v41hFdPP8y\n4LKsr1/Tr3h3fwh4qJbXEBGpmxxnxvUW835LJAetl14fkvuuB/8jJPfdbYNDclv7DwjJzUX9L+Dl\nRoOxiDSvvjIzNrO5wGqgArS5+wF5nJRIFpPeF1N8MWFQzJYu0zbMDcldX9kUkpsHr+1OiYbKY2b8\nUXdflsPriIjkK8cLePWmZQoprX+8MeaN2P5XTQzJ/c7cXUJyf8+ikNxclGiZotb3Ww48bmZTzWxM\nt88WEWmk9kq2owBqnRkf6u4LzWw48JiZvZjubPSWdJAeA3D9Dy7lC5/p6lY9kZ6Zf96kkNzr1o/o\n/kl18HJlZUhu25ZrFYqtRDPjWu8zXph+XGpm95Hs5PZkp+dobwoRidEXLuCZ2RCgn7uvTj8/Cvhu\nbmcm0o13f3+zhU9196///ruQ3D/ZjiG5v2/dNiQ3F33kAt4I4D4z63idn7r7w7mclYhIDrrZDqJQ\natmbYg6wb47nIiKSr76yZiwSqeXgT4Xk7vNvvw3JnXB3zPrny5VVIbm56CPLFCIixdZMM2Mzuxk4\nFljq7vukjw0DfgaMBOYCJ7r7X+p3miLFYUOHhOS+J2hcmd79U4qrRKXcWYo+biVje2oRkUJpb892\nFEC3M2N3f9LMRnZ6+HiSHe8haU/9BPCNHM9LpLDm3RLzJvBxYjqbLGt7MyQ3F820TNGFrtpTi4gU\nR0FmvVnUfAHP3d3MuqysUzm0NJvdDtsYkvvNX24TkntRy9YhubnoA4NxV+2p/4rKoUUkipfoAl5v\nB+OO9tSX88721CJNzwb0D8ndqiWmmmyrMt8B20xrxj1pTy0iUijNtEzR0/bUIiKF0UwzY5GiqsyZ\nGpL76mMxyxQTLeZCWislvrWtmWbGIiKlpZmxSPOasn77kNwXWlaH5M7e0OXNUsXXVp7N5bsthzaz\nm81sqZlNr3rsYjNbaGbT0uOY+p6miEgveHu2owCyzIxvBa4Fbuv0+H+7+5W5n5FIRr4mpifc3/ra\nkNzXgsqSX1pZ4u7QzbRm3MXeFCIixVeQWW8WtawZn21mnwGmAF/vagtNlUP3DRu+d27DM+c/EFME\nce/AmJ5wS4LeCZRaM82Mu3ADcAng6ccfAJ/f3BNVDt38IgZikUyafWbs7ks6PjezG4GJuZ2RiEhe\nSnQ3Ra8G445NgtIvP0nJmwFIbba68IchuSP3viYkd9ols0JyW/sPCMktNS/Pm/He7k1xmJntR7JM\nMRc4s47nKCLSO820ZtzF3hQ31eFcRHpkwKfODsm9f9X3QnJ/dM2OIbm/3WmHkNxclGgwztIDT0Sk\nnHIs+jCzr5rZDDObbmZ3mlmrmQ0zs8fM7KX0Y6/LM1UOLdJDi25+NST3N+0xc6eZaxeG5Oaiks/t\nj2a2K/AVYLS7rzOzCcBJwGiS5syXm9n5JM2Ze9UPNEs59O5m9mszm5n+VjgnfTy33wgiInWRb3fo\nFmCQmbUAg4FFJM2Zx6ffHw+c0NtTzTIzbiMp6njWzLYBpprZY8Bnyek3guSn7ZFbQnK/+c2XGp45\nqxJTBLGkbU1M7vqYrtRL1qwIyc1FxoG2ujgtNTatkQDA3Rea2ZXAfGAd8Ki7P2pmuTVnznIBbzGw\nOP18tZnNAnYl+Y1wWPq08cATaDDukyIGYpFMMq4HVxenbU76zv94YBSwArjbzE7p9BpbbM7cnR6t\nGad7VOwPTAYy/UZQOXRj9fu7Qxue+Z8PHsqGqxu/Z9Slk2LuLmixmM3lF617IyS3tWVgSG4evD23\n+4yPBF5x99cBzOxe4GB60Jy5O5kHYzPbGvg5cK67rzKzt763pd8IKodufhEDsUgm+d3aNh84yMwG\nkyxTHEGyL88acmrOnGkwNrMBJAPxHe5+b/pwbr8RRETqIqe7Kdx9spndAzxLch3tOZJJ5tbk1JzZ\nvJtyQUumwOOB5e5+btXj/wW8UXUBb5i7n7el19LMWPK0/uIvh+SOf2BYSO5v+8XsZzx747KQ3GcX\nP23dP2vL1l7zpUxjzuCzr685q1ZZZsaHAKcCfzSzaeljF5JMy3P5jSAiUhclqsDLcjfF00BXvzWO\nyPd0pIw23RpTHvzmtHUhuY9YzK1ez6z8c0jusrWrQnJz0UwbBYmIlFYzzYxFujPgsxeG5G67KiZ3\n4Y0LQnK3HTAkJHcZJZ4Z53drW93VUg6tDtEiUmyVSrajAGophwZ1iJZAf3kk5m7K7/juIbnPDInZ\nKGhiS2tIbh68mZYptlAOLSJSbM20TFGtUzk0JB2iXzCzm7vatc3MxpjZFDObMu62O2s6WRGRHslx\nP+N6q6UcOlOHaJVDS73seM0XQ3J3OiFmUrEgqBdE1G5xuSjRzLjX5dDqEC0ihddWjItzWWRpSGok\nPe9muftVVY+rQ7SE8rkxXZpfa4+5oNWfTSG5I1pL3DeiIEsQWdRSDn2yOkSLSKE10zLFFsqhH8r/\ndESya3/+jyG5+46IKYd+cmVMl+Yh/bYKyc1DU93aJiJSWs00MxYpqn57vTck9/mb5oXkzhoQc1fD\n8raYrTtzUaLBOEs5dKuZ/cHMnk/Lob+TPq7u0CJSbE1WDr0BONzd30xvcXvazH4JfAp1h5ZALcec\nEZJ7wHVfCMl9efEuIblTW8u8ZtxEM2NPdLxPGZAeTtIpdXz6+HjghLqcoYhIb7V7tqMAMpX0mFn/\n9La2pcBj7t6j7tAqhxaREO3t2Y4CyHQBz90rwH5mth1wn5nt0+n76g4tfcaIX4wLyT3roRtDcr9+\nUUyHkVwUZNabRY+K3d19BfBr4GjS7tCQVOOh7tAiUjQlWqbIUg69E7DJ3VeY2SDgY8AVwIPAaSSN\nSU8DHqjniYp0tvr0z4XkTpwSs5/xBWufC8lduWFtSO71ObyGV4qxBJFFlmWKnYHxZtafZCY9wd0n\nmtn/ou7QIlJkBZn1ZpGlHPoFkj2MOz/+BuoOLYEGnXtqSO7f/lvMTgAHbxtT5PLL118Iyc1DmW5t\nUwWeiDQvDcYi9dfyd4eH5E7p/0RI7oK2hSG5GysxW3fmojxLxjWVQ6s7tIgUmre1ZzqKoJZyaFB3\naBEpsmKMs5lkuYDnwObKoUVCbbrpkpDcvTbGvG3/v1vFNGXfbnhMZ5M8lOkCXi3l0KDu0CJSZO0Z\njwKopRxa3aEFgPbli0JyK68s7v5JdfD71piOGz9eFVP0sWRNTGeTPDTdzLhDdTm0uy9x94q7twM3\nAgfW4wRFRHqtmWbGXZVDqzu0dOg3LGaf3dZL8yiY7bmvDb8wJPeUB2L+nG9fsV9Ibh68LfoMsqul\nHPon6g4tIkXmBZn1ZlFLOXRMLaoUUtszExueeeW/T+7+SXXwq8obIbnrg6Z5aypzQnK/lceL5DgY\np9fMxgH7kExCPw/MBn4GjCSZlJ7o7r1qVtijNWORzYkYiEWy8PZsR0ZXAw+7+97AvsAsknZzk9x9\nT2BS+nWvWHIbcWPobgrJU2XO1JDcleddF5J79593C8n9Tb/VIbl3z3vAan2NpUf8Y6YxZ/ik32wx\ny8yGAtOAPbxq0DSz2cBh7r443df9CXffqzfnmnlmnN5r/JyZTUy/VndoESk0r1imo7oeIj3GdHqp\nUcDrwC3pODjOzIaQsf1cFj1ZpjiHZFreIbfpuYhIPWRdpnD3se5+QNUxttNLtQAfAG5w9/2BNXQa\n89IZc6/f/Wcq+jCz3YBPAJcBX0sfPh44LP18PPAE8I3enohIT/malSG5L8zq9eSnJjMHbQzJXdYW\n0+kjD95e80pHhwXAgqrq43tIBuMlHbf51tp+LuvM+IfAebzz2qS6Q4tIoeV1Ac/dXwNeNbOO9eAj\ngJm83X4Oamw/l6Xo41hgqbtPNbPDujhRdYeWPmPUDjHlwXu8uVNIbmtLTPl3HtxzmxkDnA3cYWYD\ngTnA50hrL/JoP5dlmeIQ4Lh0v+JWYFszu50cp+ciIvWQZ9GHu08DDtjMt3JpP5el6OMC4AKAdGb8\nH+5+ipn9F+oOLYEWnnV3SO5Va4eG5D6wOmajoA1BnT6uyuE12iu5zozrqpa2S5ej7tAiUmA5XsCr\nux4Nxu7+BMldE+oOLeF2u/2skNyjjrkrJPfprbYNyV20Nqb8Ow9NOxiLiJRJAwuMa6bBWEpr43XX\nhOS+PDDmPuOWTf1DcocFzcjzUKaZcS3l0OoOLSKF5m6ZjiLoycy4oxy6+tekukOLSGFVmu1uii7K\noUVC9Rs6KCR3nw0x+wrv1bpjSO7cfjFl53koyqw3i1rKoUHdoUWkwLzdMh1FUEs5tLpDS6i2xTH7\n7C5tibmgtaJ9Q0juiytfDcnNQ7PdTbHZcmh3P6XjCWZ2I6B2DyJSKEWZ9WZRSzm0ukNLqIEf/JuQ\n3A/PnBeSO3tlzEZBT7VXQnLzUGkvT2e5Wu4z/r66Q4tIkTXbMsVbOpVDqzu0hGo56ZyQ3Nb7vxKS\nO/r1mFnewTv0qqVbIbSX6G4KVeCJSNMq061tGoxFpGk13TKFmc0FVgMVoM3dDzCzYcDPgJEka8Yn\nuvtf6nOaIsUxYf6uIbkP94/ZPe2F1TEXLPNQpmWKnixCfdTd93P3jp3u1R1aRAqt0t4v01EEtSxT\nqDu0hLLWISG5Z94f00eh9bgJMbnblnc1s0SrFJlnxg48bmZTzWxM+pi6Q4tIobW7ZTqKIOuvvEPd\nfaGZDQceM7MXq7+p7tASYdFRY7p/Uh08vvRdIbkLW2L+95m3obyXgsp0N0WmmbG7L0w/LgXuAw4k\n7Q4NoO7QIlJE7RmPIsiyUdAQoJ+7r04/Pwr4LvAg6g4tqbanGr+euX7tgIZnAoxL5iYN9+KymA17\nNpW4HNopz8w4yzLFCOA+M+t4/k/d/WEzewZ1hxZiBmKRLNpKtEyRZaOgOcC+m3lc3aEFgJYPn8jy\nT//V7ql1d8f6mPt99xmwMSR3x2Exd4/MXvdaSG4emm1mLLJFEQOxSBZFWQ/OQoOxiDStppsZd1EO\nfTFwBvB6+rQL3f2hepykFNuwe24Oyf1WSCps+N65IbmP/zSmw8j4QeUZ0Dpr1pnxR919WafH1B1a\nRAqr0mwzY5EtaZv5ZEju50+9PyR3xobXu39SHVR8SUjum2vWheTmoURdl2oqh4YM3aFFRKK0Y5mO\nIqilHDpTd+h08B4DcP0PLuULnzk5lxOX4thww/iQ3GM3xdzaNrWyICR3zsrF3T9J3qFM+y9kGoyr\ny6HN7D7gQHd/673plrpDa28KEYnSVBfwuiqHVndo6TDkmptCcj950yUhuZUf7R2S+9Quo0JyH109\nOyQ3D+1WjCWILGoph/6JukOLSJGVaVeNWsqh1R1aRAot77spzKw/MAVY6O7H5tl+Tre2SWmteXxO\nSO7/a9kuJHdZ29qQ3FUbYnLzUIc7Jc4BZgEdFTgd7ecuN7Pz06971fGoGM2fRETqwDMeWZjZbsAn\ngHFVDx9P0naO9OMJvT3XrOXQ26UnsA/JuX8emI26QwtQmR9z7XbGC8NDcltaN4XkTlvxSkhumWVd\npqi+BTc1Nr0TrNoPgfOAbaoey9R+LousM+OrgYfdfW+S9eNZqDu0iBRc1k4f7j7W3Q+oOt4xEJvZ\nscBSd5/aVZa792Si/Vey3No2FPgI8Nk0cCOw0czUHVoA2HTzj0Nyh28dc9vSm5ticnffeqeQ3KXr\nV4Tk5qGS33+qQ4DjzOwYoBXY1sxuJ20/5+6La20/l2VmPIpkZ7ZbzOw5MxuX3m+c2/RcRKQe8uqB\n5+4XuPtu7j4SOAn4lbufwtvt56DG9nNZ1oxbgA8AZ7v7ZDO7mk5LElvqDq1y6ObXevG1+Po1Dc8d\n9NTXGp4JsN2Kbbp/Uh08vnJRSG6ZNaAC73Jyaj+XZTBeACxw98np1/eQDMaZpucqh25+EQOxSBb1\naIHn7k+QLMvm2n4uS9HHa2b2qpnt5e6z0+CZ6aHu0IK1xvRmGzFmr5Dcfa+IWUO9s9IWkltmTbU3\nReps4A4zGwjMAT5Hst6s7tAiUlhNVQ4N4O7TgAM28y11hxaRwirT5vIqh5aatf3u3pDc+y6LWS5Y\nNDBmueDDw0eH5P5x1byQ3Dw04zKFiEjpaDCWPqXl4E+F5O5QmRKSuzpoJbLMM9QoZbp9K1M5tJlt\nZ2b3mNmLZjbLzP7BzC42s4VmNi09jqn3yYqI9ES7ZTuKIOvMuGNvik+nd1QMBv4J+G93v7JuZyey\nBYfdf3xI7ntPvzUk98y2mDXjC219SG4emupuii3sTVHfMxMRqVF7iRYqssyMq/em2BeYSrLBMsDZ\nZvYZkp3vv765LTRVDi310n/Pvw/J3f3umK07h5x5aUjuiFe3D8nNQ5ku4GVZM+7Ym+IGd98fWENS\nDn0DsAewH7AY+MHm/uXqrek0EItII+W5uXy9ZRmMN7c3xQfcfYm7V9y9HbgROLBeJyki0ht57drW\nCL3em6Jjk6D0aZ8EYto9SJ9VmdPlPt919ZevXR+Se+2CnUNyV/iqkNw8tG1+M8lCqmVvih+Z2X4k\ns/y5wJl1OUMRkV4qz1Bc294Up+Z/OiLZ9d/jgyG52533iZDcbzzwSEju3Y+Ut29EUZYgslAFnog0\nrWa7tU2kkNoXvRSSO/mMyd0/qQ6e3ypmhvr7/m+G5J6ew2uUZyjOcDeFme1VVfI8zcxWmdm5ZjbM\nzB4zs5fSj+W9GVFEmlKz3U0xm+ReYsysP7AQuI/kXuNJ7n65mZ2ffq3u0NI4Q4aGxPYPukJ/3boX\nQ3KXbyjv3RSVEs2NM20UVOUI4M/uPg84HhifPj4eOCHPExMRqVVTzYw7OQm4M/18RNV9xq8Bm13Q\nUjm01Ev7n58NyR31nuUhuV9ZtHdI7i9a3wjJzYOXaGaceTBO7zE+Drig8/fc3c02/95N3aFFJEpR\nZr1Z9GSZ4uPAs+6+JP16iZntDJB+XJr3yYmI1KIdz3QUQU+WKU7m7SUKgAeB04DL048P5HheIt1b\nsSwkduj+A0Nyd5sXM89rHVjeO2CLMcxmk+lP2cyGAB/jnSXPlwMTzOx0YB5wYv6nJyLSe20lGo6z\nlkOvAXbo9NgbJHdXiMTYaZeQ2C/9cnBI7oqBMbeY/fK150Jy89CUF/BERMqmTBfwNBhLac3/4s9D\nct8dVGw6wHpaFpCP92xb3o2CmmpmbGZ7AT+remgP4CJgO+AMkpZMABe6+0O5n6GISC811cx4C+XQ\nn0PdoSXQHk9fF5L7zVu/F5K76OaY4ovVvntIbh4q3kQz407eKodWd2gRKbqi3EOcRS3l0KDu0JJq\nX74oJPdbRzV+dnzbimkNzwTY0LYpJHfNpvUhuXko05px5isCVeXQd6cPqTu0AH1rIJZyadaNgt5R\nDl1VFo2Z3QhMzPncRERq0qzLFO8oh1Z3aOnw+qkXheR+aENMt+SnhrwrJHdNJWa5oOJFmTv2XF7L\nFGa2O3Abye6UDox196vNbBjJ3WYjSRozn7i55dosMi1TVJVD31v18PfN7I9m9gLwUeCrvTkBEZF6\nqbhnOjJoI7kuNho4CDjLzEbzdpONPYFJ6de9Uks5tLpDCwAjfjEuJDeqm8ExP455J7Dq4YUhudcu\niHkHkoe8linSVYDF6eerzWwWsCtJk43D0qeNB56glx2PYkp6REQaIOsFPDMbY2ZTqo4xXb2mmY0E\n9gcmk7HJRhYqh5bS+u37YlouXjEwplvyzLUxRR8LV/8pJPeyHF4j65pxdROMLTGzrYGfA+e6+6rq\neostNdnIIuua8VfNbIaZTTezO82sVd2hRaTo8txc3swGkAzEd7h7x/Wz3JpsZNmbYlfgK8Bod19n\nZhNIij9Go+7QEuiQGVeE5H7/Q+eE5P5h0PtDcudsU+K7KXIqh7ZkCnwTMMvdr6r6Vm5NNrKuGbcA\ng8ysBRgMLELdoUWk4Cp4piODQ4BTgcPNbFp6HEMyCH/MzF4Cjky/7pUsGwUtNLMrgfnAOuBRd3/U\nzNQdWkQKLce7KZ4GutqQJ5cmG1mWKbYnmQWPAlYAd5vZKdXPUXfovs3XrwnJfeXI80Jy7/KYyyNz\n+sdcOHxy5UshublcwGuyXduOBF5x99cBzOxe4GDShWt3X6zu0CJSRM1WDj0fOMjMBpMsUxxBskvb\nGtQdWgBrHRKSu8uRMdu4/vO9MWXJl7W0heSWWZl2bcuyZjzZzO4BniUpCXyOZNlha9QdWkQKrOk2\nl3f3bwPf7vTwBtQdWgK1XnxtSO6Ip84MyX3/yp1Cch9fH7NWnYdmW6YQESklDcYiTWzK0uEhuc+2\nLA/J3XPoLiG5eSjT3RS1lENfbGYLO90ALSJSGHmWQ9dbLeXQoO7QEmjZCaeH5P7vVr3emKsmTyyZ\nGZJbaa+E5Oahqe6mqHreIDPbxNvl0CPrdVIiInkoU5eSbpcp3H0h0FEOvRhY6e6Ppt8+28xeMLOb\nu9q1rXqf0HG33bm5p4iI1IW7ZzqKoJZy6BuAS0j6QV1C0h36853/fZVDS70MHhXTG2GX+THXvftb\nUC+IEregKMp6cBZZ/pjfKod2900kffAOdvcl7l5x93bgRuDAep6oiEhPecZ/iqDX5dDqDi0d2h65\nJST3B0/FdGle0C+mHHrrga0hucvWrgrJzUN7QZYgsqilHHqcme1HskwxF4gpSxIR6UJRZr1ZWCMX\nr7VmLM3gifddEJL73ZZlIbnL22LKoWcsmVzzTlB7D/9QpjHnxaXPxOw6VUUVeCLStJpqmUKkqCa8\n/6KQ3KmDYiZRs5a/GpJbZmVapshaDn1OWgo9w8zOTR9Td2gRKbR290xHEXQ7GJvZPsAZJLeu7Qsc\na2bvJekGPcnd9wQmpV+LiBRGs93a9n+Aye6+FsDMfgN8iqQQ5LD0OeOBJ4Bv5H+KUnQvBrWuX9wy\nNCT3VV8Xkjtsq21DcstUUtxZxcuzr0aWZYrpwIfNbIf0XuNjgN2BzN2hVQ4tIhGaqhza3WeZ2RXA\noyR976YBlU7PUXfoPmzvZ64OyX3vQzeG5M6/aE5I7msDY94JzFpd3guHzVYOjbvf5O4fdPePAH8B\n/kTaHRpA3aFFpIiaamYMYGbD3X2pmb2bZL34IJKNg9QdWsLM/960kNxX2vuH5EYVXwwduHVIbh6K\ncqdEFlnvM/65me0AbALOcvcVZnY56g4tIgVWlDslssjaHfrDm3nsDdQdWlJtQeu3Ef6lsl1I7u0D\n20Jy/7yuvCuQZboTRBV4UrO+NBBLuRRlPTgLDcZSs5ZjzgjJfffIJ0Nyh199a0juoim7heQOGTwg\nJDcPZVozrqUcWt2hRaTQmupuik7l0BuBh81sYvptdYcWkcIq033GtZRDiwDQ9uzDIbnr/2dCSO7/\nBC0XTGV1SO5rlZhb6vJQlFlvFrWUQ4O6Q4tIgVW8PdNRBJk6faT3En+JpBx6BrAB+E9gGW93h97Z\n3f+qO3Q1lUNLnt48c4t/3erm8d/tGpJ7Vb9FIbltQZvtPLPoyZo3jh406D2Zxpx16+aFd/rodTm0\nukOLSNE11QU82Hw5tLpDS7SWnWK6JUfNUKctj9mgaFMlptgkD3lW4JnZ0cDVQH9gnLtfntuLU1s5\n9DXqDi0iRZbXrNfM+gPXAR8DFgDPmNmD7j4zlwBqK4c+Na+TEBGphxyLPg4EXnb3OQBmdhdJg43c\nBuPMayrRBzCmr+T2pZ+1r+X2pZ81Mrc35wlMqTrGdPr+p0mWJjq+PhW4Ns9zyHQBryDG9KHcvvSz\n9rXcvvSzRub2iLuPdfcDqo6xjT6HMg3GIiJRFvJ2fQXAbuljudFgLCLSvWeAPc1slJkNBE4CHswz\noEy7tjX8bUNgbl/6Wftabl/6WSNzc+XubWb2ZeARklvbbnb3GXlmZKrAExGR+tIyhYhIAWgwFhEp\ngEIOxukucEvNbHrVY8PM7DEzeyn9uNld4mrI3N3Mfm1mM9NN9M9pUG6rmf3BzJ5Pc7/TiNw0o7+Z\nPdexP3WDMuea2R/ThgRTGpi7nZndY2YvmtksM/uHBvy33auq+cI0M1tlZuc2IPer6d+l6WZ2Z/p3\nrBF/xptrQlH33GZRyMEYuBU4utNj5wOT3H1PYFL6dZ7agK+7+2jgIOAsMxvdgNwNwOHuvi+wH3C0\nmR3UgFxEUnZ/AAADQklEQVSAc4BZVV83IhPgo+6+n7sf0MDcq4GH3X1vYF+Sn7uuue4+O/059wM+\nCKwF7qtnrpntCnwFOMDd9yG52HRSPTPT3OomFPsCx5rZe+ud21SiK1+2UBEzEphe9fVskm06AXYG\nZtc5/wGSOvSG5QKDgWeBv693Lsl9kpOAw4GJjfozJtnHZMdOj9X7Zx0KvEJ6wTri7xRwFPDbeucC\nuwKvAsNI7paamGbX+8/4X4Cbqr7+FnBeo/+/LfNR1Jnx5ozwt3eJew0YUa8gMxsJ7A9MbkRuulww\nDVgKPObujcj9Icn/LNU7azfiz9iBx81sqpl1VGfVO3cU8DpwS7osM87MhjQgt9pJQEd3hbrluvtC\n4EpgPrAYWOnuj9YzM9VVE4pG/hmXWpkG47d48mu2LvfkmdnWwM+Bc919VSNyPdkXej+S2eqB6Vu+\nuuWa2bHAUnefuoVzqtef8aHpz/pxkqWgjzQgtwX4AHCDu+9P0iThHW+X6/x3aiBwHHB35+/V4b/t\n9iQb2IwCdgGGmNkp9cxMX3MWcAXwKPAwMA2odHpO3f6Mm0GZBuMlZrYzQPpxad4BZjaAZCC+w93v\nbVRuB3dfAfyaZL28nrmHAMeZ2VzgLuBwM7u9zpnAWzM33H0pyfrpgQ3IXQAsSN9xANxDMjg36r/t\nx4Fn3X1J+nU9c48EXnH31919E3AvcHCdM4HNN6FoRG6zKNNg/CBwWvr5aSRrurkxMwNuAma5+1UN\nzN3JzLZLPx9Esk79Yj1z3f0Cd9/N3UeSvH3+lbufUs9MADMbYmbbdHxOspY5vd657v4a8KqZ7ZU+\ndATJ1od1za1yMm8vUVDn3PkkzR8Gp3+njyC5WFn3n9XMhqcfO5pQ/LQRuU0jetF6cwfJX9zFJJvZ\nLwBOB3YgueD0EvA4MCznzENJ3kK9QPIWaxrJule9c98PPJfmTgcuSh+va25V/mG8fQGv3j/rHsDz\n6TED+GajflaSO1WmpH/O9wPbNyh3CPAGMLTqsXr/OX+H5Bf6dOAnwFYN+lmfIvkl9zxwRCP/HjfD\noXJoEZECKNMyhYhI09JgLCJSABqMRUQKQIOxiEgBaDAWESkADcYiIgWgwVhEpAD+P+H1kQY5NUQf\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2aab0a2e6b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(cnvs.T, xticklabels=quality, vmax=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, 100 or less CNVs might not be a bad value to start with. Now we need to split which ones are de novo and inherited."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO:\n",
    "\n",
    "* split CNVs into denovo and inherited\n",
    "* better to run PCA within groups? are controls biasing the results? check literature. might need ot run other tools instead"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
